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ABSTRACT 

We present a numerical investigation of dead, or relic, radio galaxies and the environ- 
mental impact that radio galaxy activity has on the host galaxy or galaxy cluster. We 
perform axisymmetric hydrodynamical calculations of light, supersonic, back-to-back 
jets propagating in a /3-model galaxy/cluster atmosphere. We then shut down the jet 
activity and let the resulting structure evolve passively. The dead source undergoes 
an initial phase of pressure driven expansion until it achieves pressure equilibrium 
with its surroundings. Thereafter, buoyancy forces drive the evolution and lead to the 
formation of two oppositely directed plumes that float high into the galaxy/cluster at- 
mosphere. These plumes entrain a significant amount of low entropy material from the 
galaxy/cluster core and lift it high into the atmosphere. An important result is that 
a large fraction (at least half) of the energy injected by the jet activity is thermalized 
in the ISM/ICM core. The whole ISM/ICM atmosphere inflates in order to regain 
hydrostatic equilibrium. This inflation is mediated by an approximately spherical dis- 
turbance which propagates into the atmosphere at the sound speed. The fact that such 
a large fraction of the injected energy is thermalized suggests that radio galaxies may 
have an important role in the overall energy budget of rich ISM/ICM atmospheres. 
In particular, they may act as a strong and highly time-dependent source of negative 
feedback for galaxy/cluster cooling flows. 

Key words: cooling flows — galaxies: jets — hydrodynamics — radio galaxies — 
X-rays: galaxies: clusters 



1 INTRODUCTION 

There is now overwhelming observational and theoretical 
evidence that powerful radio galaxies possess highly colli- 
mated and relativistic twin jets of matter that emerge from 
a central active galactic nucleus (AGN). In the powerful 
Faranoff-Riley type 11 radio galaxies (FR-II; Fanaroff & Ri- 
ley 1974), the jets are thought to remain relativistic along 
most of their length before terminating in a series of shocks 
resulting from interaction with the surrounding material. A 
substantial amount of theoretical and numerical (hydrody- 
namic) work suggests that, after passing through the ter- 
minal shock, the spent jet material inflates a broad cocoon 
which encases the jets (e.g., see Scheuer 1974; Norman et 
al. 1982; Lind et al. 1989; Begelman & Ciofii 1989; Ciofli & 
Blondin 1992; Hooda, Mangalam & Wiita 1994). 

In the early life of a powerful source, the cocoon is prob- 
ably overpressured (Begelman & Ciofli 1989) with respect 
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to the surrounding ambient material (either the interstel- 
lar medium [ISM] of the host galaxy, or the intracluster 
medium [ICM] of the host cluster). The cocoon then un- 
dergoes supersonic pressure-driven expansion into the am- 
bient medium, sweeping the ambient medium into a shocked 
shell. A contact discontinuity separates the relativistic co- 
coon material from the shocked ambient material. From the 
onset of activity, Kelvin-Helmholtz (KH) instabilities act to 
shred the contact discontinuity. If the ambient medium has 
a density profile characterized by p oc r^" where a > 2, 
the contact discontinuity will accelerate and hence will be 
Rayleigh- Taylor (RT) unstable. However, in the more phys- 
ical case of Of < 2 the contact discontinuity will possess an 
initial deceleration, which is sufficient to stabilize the con- 
tact discontinuity against RT modes. Only at later times, 
once the cocoon comes into approximate pressure balance 
with its surroundings, will the deceleration of the contact 
discontinuity no longer exceed the local gravitational accel- 
eration, thereby rendering it RT unstable. Eventually, AGN 
activity will cease and the jets will turn off. Although this 
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stage of a radio galaxy's life has been little studied, it seems 
likely that the combined action of the KH and RT instabili- 
ties will transform the remnant of the cocoon into 'plumes' 
that rise in the potential of the galaxy/cluster under the ac- 
tion of buoyancy forces (Gull & Northover 1973; Churazov 
et al. 2000; Briiggen & Kaiser 2000). 

Observationally, radio galaxies arc often seen to bo as- 
sociated with galaxy or cluster cooling flows. This raises 
an obvious and interesting set of questions. To what extent 
is radio galaxy activity a natural response to the cooling 
flow phenomenon? Do radio galaxies act back on their en- 
vironment to partially offset the cooling flow, i.e., are radio 
galaxies a dramatic manifestation of Ic Chatelier's princi- 
ple? Is radio galaxy activity a crucial component for our un- 
derstanding of galax;y clusters? Since it takes only a small 
fraction of the cooling flow mass to fuel a powerful AGN, 
it is easy to sec how radio galaxy activity can result from a 
cooling flow. However, the impact of a radio galaxy on its 
environment is much less well understood. 

In this paper, we use hydrodynamic simulations to 
study the evolution of a radio galaxy situated at the center 
of a galaxy/cluster. In particular, wc follow the evolution of 
the source at times after the jets have terminated. We ex- 
amine the interaction of the radio galaxy with the ISM/ICM 
and assess the environmental impact of the AGN activity. 
Section 2 presents our set of hydrodynamic simulations and 
focuses on some numerical issues. In Section 3, wc briefly dis- 
cuss the 'active' phase of the source. Section 4 presents our 
results for dead sources (i.e., the passive phase after the jets 
have turned off), and Section 5 discusses some astrophysical 
implications as well as the limitations of our calculations. 
Section 6 draws our conclusions together. 



2 THE HYDRODYNAMIC SIMULATIONS 

2.1 Basics and simulation setup 

We model a source with back-to-back supersonic jets sit- 
uated at the center of a spherical gaseous atmosphere. To 
permit parameter studies of various high resolution simula- 
tions, we make the assumption of axisymmetry and neglect 
magnetic flclds. Wo discuss these assumptions in Section 5. 
The simulations were performed using the ZEUS-3D code 
(Stone & Norman 1992a, 1992b; Clarke, Norman & Fiedler 
1994) in its 2D (axisymmetric) pure hydrodynamic mode. 
ZEUS-3D is an explicit Eulerian finite difference code, for- 
mally of second order, which uses an artificial viscosity to 
reproduce shocks. Major advantages of this code are its fiex- 
ibility in the choice of the computational grid, and the fact 
that it has been extensively tested. The simulations were 
performed on a Sun Microsystems Ultra 60 workstation. 

Unlike many previous simulations, wc model both jets 
and make no assumption regarding reflection symmetry in 
the plane normal to the jets. This allows us to study phe- 
nomena associated with the collision of the two backflows. 
Our simulations were performed in spherical polar coordi- 
nates {r, 0, (f>), and the computational domain was the region 
r £ (0.1, 10), except for one run (Run 4 below) which had 
the domain r £ (0.1, 30). All gas within the calculation was 
assumed to have an adiabatic index of 7 = 5/3. The ambient 
galax;y /cluster gas was assumed to be initially in an isother- 



mal configuration with a (adiabatic) sound speed Cs = 1 and 
a density profile given by a /3-model with /3 = 0.5: 

^^"■^ " [1 + (r/ro)2]3/4 ' 

where the core radius was set to ro = 2. The gravitational 
potential, <I>, was set so as to keep this ambient material in 
hydrostatic equilibrium 

— = 1^ (2) 

dr p dr ' 

which gives (to within an additive constant) 

*=fln(p). (3) 

It is assumed that the gravitational potential is dominated 
by a stationary distribution of background dark matter, i.e., 
the self-gravity of the gas is negligible. This reasonable as- 
sumption prevents us from having to solve Poisson's equa- 
tion as part of the hydrodynamic simulations. 

We inject back-to-back jets into this background 
medium via the use of inflow boundary conditions at r = 0.1. 
The jets are initially conical with a 15° half opening angle, 
and in pressure equilibrium with the ambient medium. They 
are given an initial density of pjet = 0.01 and Mach number 
(with respect to the sound speed of the injected jet material) 
of TVf = 10. As discussed in the Appendix, these parameters 
are chosen so that the Kelvin-Helmzholtz growth rate at 
the cocoon-ICM contact discontinuity in this non-relativistic 
simulation is approximately the same as the corresponding 
growth rate in the real relativistic situation. At time t = 1, 
the jet activity is stopped and the inner boundary is made 
reflecting. The source is then allowed to evolve passively. 

We can relate quantities within the code to physi- 
cal quantities by fixing the parameters of the background 
medium. Suppose we set ro = 100 kpc, cism = 1000 km s"'^ 
and a central number density of no = 0.01 cm^"* — values 
representative of a rich galax;y cluster. Then, one code unit 
of time corresponds to 50 Myr. The total kinetic luminosity 
of the source is then 9.3 x lO^^ergs"^ . These parameters 
are relevant for powerful sources in rich galaxy clusters such 
as Cygnus A (Arnaud et al. 1984; Reynolds & Fabian 1996) 
and 3C 295. 

Unless otherwise stated, this paper shall use the above 

scalings when converting our simulation results into quanti- 
ties that can be compared directly with real systems. How- 
ever, in the present universe at least, many radio galax- 
ies are found in smaller galaxy clusters or poor groups. 
Taking objects like Hydra A and Virgo A as prototypes, 
the typical relevant ISM/ICM parameters are ro — 10 kpc, 
Cism = 500 km s^^ and no = 0.1 cm"'' (e-g-, see the detailed 
analysis of the Chandra X-ray Observatory data for Hydra A 
by David et al. 2000). In this case, one time unit corresponds 
to 10 Myr and the total kinetic luminosity of the source is 
1.16 X 10'*'* erg s^^ . Of course, it must be appreciated that 
these /3-models are very crude parameterizations of the real 
ISM/ICM density structure. 

2.2 The need for adequate numerical resolution 

During the active phase, the large scale evolution and prop- 
erties of the cocoon are strongly infiuenced by much smaller 
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scale shock structures that appear in the jet. As wiU be dis- 
cussed in more detail below, a conical shock appears in the 
jet which sprays the jet thrust over a large working surface 
(which expands almost self-similarly) . Failure to resolve this 
shock structure in a simulation leads to the jet thrust be- 
ing deposited over a very small area which leads to the jet 
rapidly 'drilling' a narrow cavity through the ICM. 

Initially, we performed simulations with rout ~ 10 at 
three different resolutions: iir x ne — 300 x 300 (Run 1), 
600 X 600 (Run 2), and 1200 x 1200 (Run 3). The resolu- 
tion was slightly enhanced at small r and along the jet axes. 
We find that Run 1 does indeed fail to resolve the jet-shock 
structure. On the other hand, Run 2 and Run 3 both well 
resolve the internal shocks within the jet and produce co- 
coons that expand in almost a self-similar manner, at least 
during the early active phase. The similarity of these two 
simulations gives us some confidence that we have achieved 
a degree of numerical convergence. 

We note one additional numerical technicality. Our sim- 
ulations possess a relatively high dynamical range in r. How- 
ever, to maintain adequate resolution on large scales, we 
choose not to use a logarithmic ("scaled") grid for the r 
variable. Instead, we space our grid using a more general ge- 
ometric progression (i.e. "ratioed" scaling). While allowing 
us to maintain adequate resolution on the scales that matter 
to our problem, it has the undesirable effect of forcing highly 
elongated grid elements to exist in the very central regions 
of the simulation. One might then worry about the effects 
of an effective anisotropic artificial viscosity in this inner re- 
gion. However, this is acceptable since we are not interested 
primarily in the detailed hydrodynamic properties of these 
inner regions. One can view the innermost region of our grid 
as a buffer zone, which connects smoothly to the interesting 
region of the grid, while protecting it from the hard inner 
boundary. Experiments with high resolution simulations us- 
ing a logarithmic r-grid give us confidence that this is not a 
problem. 

For the rest of this paper, we discuss results from Run 2 
(i.e., Tout = 10 with a grid size of rir x ng = 600 x 600), and 
one other simulation. This final simulation (Run 4) has a 
larger outer radius (rout = 30) and a larger computational 
grid (rir xng — 1200 x 600). In the spatial region covered by 
both Run 2 and Run 4, the resolutions are approximately 
(but not exactly) the same. 



3 THE ACTIVE PHASE 

There is already a substantial body of numerical work on 
the hydrodynamics, and magnetohydrodynamics, of active 
radio galaxy jets. Thus, we shall discuss only briefly our 
simulation results for the active phase of our source. In a 
companion paper (Reynolds, Heinz & Begelman 2000) we 
focus on the X-ray appearance of such active sources. 

Our results during the active phase are very similar to 
those obtained by Lind et al. (1989), who performed ax- 
isymmetric hydrodynamic simulations of a single jet prop- 
agating through a uniform medium. The initially conical 
jet is focused by the external pressure and suffers a series 
of oblique shocks. Near the working surface with the back- 
ground medium, a strong annular shock deflects most of 
the flow into a wide fan. Thus, the jet thrust is deposited 



over a large area, thereby reducing the advance speed of 
the jet below the value that would be obtained if the jet 
just 'drilled' into the background medium. This is the 2- 
dimensional equivalent of the so-called "dentist drill effect" 
(Scheuer 1974). After passing through the annular shock, the 
spent jet material enters a mildly supersonic backflow which 
gradually decelerates through a series of weak shocks. Power- 
ful Kelvin-Helmholtz (KH) instabilities exist along the con- 
tact discontinuity between the backflowing jet material and 
the shocked/compressed ambient material. These KH insta- 
bilities, which become more powerful as the backflow decel- 
erates toward the region where the two backflows collide, act 
to mix ambient material into the cocoon. These instabilities 
can be seen clearly in Fig. 0, which shows specific entropy 
maps {a = p/p^^^; note that we neglect logarithms and con- 
stants in our definition of the specific entropy) for Run 2 at 
times t — 0.5, t — 1.0, t = 3.0 and t = 5.0. During the active 
phase, the structure expands approximately self-similarly, at 
least as long as the cocoon is still sufficiently overpressured 
as to undergo supersonic lateral expansion into the ambient 
material. 

It is worth noting explicitly that the backflow in the 
cocoon has a speed comparable to the sound speed of the 
cocoon material. In code units, the advance speed of the 
working surface of the jet is ~ 3 (i.e., a Mach number of 3 
with respect to the sound speed of the undisturbed ambient 
material). This should be compared with a backflow speed 
of ~ 20 — 30 which exists in a large portion of the cocoon. 
Translated into terms applicable to real radio galaxies, the 
backflow speed would be a significant fraction of the speed 
of light (since the spent jet material is a relativistic plasma 
with a sound speed of c/\/3). 



4 THE PASSIVE PHASE 

As discussed in Section 2, we turn off the jets after t = 1 
and allow the existing structure to evolve. Both Run 2 and 
Run 4 display very similar behavior during the passive phase 
although, of course, Run 4 can be followed to much later 
times due to the larger spatial domain of the calculation 
(see Fig. 2). We will describe our findings based on Run 4. 
It should be noted that Briiggen & Kaiser (2000) and Chu- 
razov et al. (2000) have also simulated the late evolution of 
radio galaxies. Their simulations make the assumption that 
the source evolution before the jets turn off is governed by 
simple supersonic expansion of the lobes and consequently, 
their initial conditions for the passive phase are based on rel- 
atively simple geometries; namely, spherical and ellipsoidal 
bubbles placed in isothermal atmospheres. However, turbu- 
lence and powerful circulatory motion seems to be present 
at all times inside the active cocoon, which leads to signif- 
icant excitation of the Kelvin-Helmholtz instability before 
the source turns off. Whether, and by how much, the co- 
coon is overpressured at the moment that the jet activity 
ceases can be estimated by the simple analytic models of 
Begelman & Ciofh (1989). 

Thus, we extend the work of Churazov et al. (2000) and 
Briiggen & Kaiser (2000) by simulating the active phase of 
the source and then shutting the jets off, rather than mod- 
elling the late stages of evolution by letting a static bubble 
evolve under the action of buoyancy. We also simulate a 
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Figure 1. Specific entropy maps for t = 0.50, t = 1.0, t = 3.0 and t = 5.0 for Run 2. As the same greyscale levels are used for all four 
panels, one can clearly see the decrease in the specific entropy of the cocoon it rises and entrains ambient thermal plasma. Note the 
change of scale in the first three of these plots. 



substantially larger portion of the cluster atmosphere while 
maintaining the required high resolution, thereby aUowing 
us to track the source evolution at later times. We consider 
these aspects to be important extensions of these previous 
works. 



4.1 The stages of an inactive source's life 

As would be expected from simple qualitative arguments, 
turning off the jets results in the almost immediate disap- 
pearance of internal shocks in the jet channels, which sub- 
sequently collapse. The very rapid backflow motions in the 
cocoon slow down accordingly, due to the lack of high- inertia 
material and the disappearance of the terminal shock, which 
previously re-directed the flow. While the radial lobe expan- 
sion during the active phase was governed mostly by the jet 
thrust, the absence of ram pressure from the jets in the in- 
active phase results in a rapid slowing of the cocoon head to 
subsonic speeds. If the cocoon is still signiflcantly overpres- 
sured when the jet activity ceases, this subsonic stage will 
be preceded by a Sedov expansion phase. 

During the early inactive phase of our simulations, the 
source is still slightly overpressured with respect to the am- 
bient medium and we do indeed find a Sedov phase, during 
which the cocoon expands in order to achieve pressure bal- 
ance. When pressure balance is reached (at t ~ 1-5), the 
driving force on the shocked/compressed ICM shell ceases 



and this shell rapidly decelerates to subsonic speeds. Once 
the lobes reach rough pressure equilibrium with their sur- 
roundings, the dynamics of the lobes closely resemble a ris- 
ing bubble. In the initial stages of this late inactive, or buoy- 
ant, stage, the denser IGM begins to settle back into the 
center, squirting the light cocoon material along the major 
axis of the cocoon. Eventually, the ICM core will reform and 
the cocoon will be pinched off into two separate, buoyantly 
rising plumes. 

Material that had been entrained previously is now car- 
ried along by each plume, partly due to its own inertia and 
partly due to the significant circulation inside the plume. 
This vorticity continues to drag ambient material into the 
bubble at the plume's back end, where strong eddies pro- 
duce a jelly-fish- like structure (Fig. At the head of the 
plume, a dense cap of swept up material forms, which is 
pushed ahead. 

Continuing with the bubble analogy, we can write the 
buoyant velocity Vc of the plume as 



4rbubg(pa — /9c) 



(4) 



where rbub the radius of the bubble, g is the local gravita- 
tional acceleration, and pa, and pc are the densities of the am- 
bient and bubble materials respectively. The top left panel 
of Fig. ^ shows the volume averaged plume height (i.e., its 
distance from the center). The plumes have approximately 
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Figure 2. Specific entropy maps for t = 5.0, t = 10.0, t = 15.0 and t = 20.0 for Run 4. It can be seen that low entropy material, from the 
center of the galaxy/cluster atmosphere, is being lifted by the buoyant plumes. Note the change in scale between all of the four panels. 




Figure 3. Entropy map of a buoyant plume. The arrows corre- 
spond to the 2D velocity field. 



constant velocity of Vc ~ 0.8. Together with the cocoon size 
of rhnh{t = 4) ~ 1 and the gravitational acceleration at that 
position, this indicates that C is sUghtly less than unity. 



4.2 Energetics and Thermodynamics 

One goal of these simulations is to determine a set of phys- 
ical parameters that are only weakly dependent on the de- 
tails of the simulations (most importantly, they should be 
weakly dependent on the numerical resolution and the di- 
mensionality). In particular, we are interested in the over- 
all matter distribution, energetics, and thermodynamics of 
the system. For this discussion, it is interesting to sepa- 



rate the system into different structural components, i.e., 
the cocoon/plumes, and the ambient medium (including the 
shocked/compressed shell that bounds the radio source). 

Defining these regions is not trivial. Since ZEUS-3D is 
an Eulerian code we cannot, a priori, separate cocoon plasma 
from ambient plasma by tracking fluid elements. Due to 
the turbulent cocoon surface, a simple geometric identifi- 
cation of the cocoon is also not an option. However, since 
the jet material is of very high entropy compared to the 
background gas (especially once it has passed through the 
thermal shock), we use the specific entropy s = p/p^^^ as a 
discriminant. 



4-2.1 The entropy threshold method and numerical mixing 

For the simulation parameters we chose, the background ma- 
terial has specific entropy (in code units) 



P 

„5/3 



= |[l + (r/ro)1^/^ 



(5) 



The injected material, on the other hand, has entropy 
p/p^/^ ~ 1300 and cocoon material, which has gone through 
the terminal shock, can have entropy as high as s ~ 10^. We 
chose to define the cocoon as the volume inside of which 
the entropy satisfies Sc > 10. At all but the very earliest 
times in our simulations, the shock driven into the ambient 
material is too weak to raise its entropy above this thresh- 
old. Conversely, an examination of entropy maps from the 
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Figure 4. Simulations diagnostics as a function of time. Top left panel : Volume averaged height of the plume as a function of time. 
Top right panel : Ratio of mass entrained in the cocoon/plume against mass injected by the radio jets. Bottom left panel : Energetics 
of the ambient ICM during the evolution of the radio galaxy. Shown here are the internal energy (ABint,ambi solid line), kinetic 
energy (AEjjin, amb) dotted line), and gravitational potential energy (Aii/pot^ambi dashed line). Bottom right panel : Energetics of the 
cocoon/plumes during the evolution of the radio galaxy. Shown here are the internal energy (A_Eijit,coci solid line), kinetic energy 
(Ai?]jin coci dotted line), and gravitational potential energy (Ai?pot,coc, dashed line). 



simulations shows that, until about t ~ 8, the contact dis- 
continuity between the cocoon and ambient material is sharp 
and well tracked by the Sc = 10 contour. In fact, until t ~ 5, 
the location of the contact discontinuity is insensitive to the 
chosen entropy threshold, provided that threshold is in the 
range 10 < Sc < 100. For times t > 8 numerical mixing 
reduces substantial parts of the cocoon/plumes below the 
Sc = 10 threshold. At these late times, identification of a 
well-defined contact discontinuity separating the plume from 
the ambient material becomes problematic. 

To better illustrate the effects of numerical mixing at 
later times, we calculate the entropy distribution of Run 4. 
Figure ^ shows this function at different times. Due to nu- 
merical mixing the second peak at high entropy (i.e., the 
cocoon material) diffuses to lower entropies. It is clear that 
after a certain amount of time an entropy cut is not a reason- 
able estimator for cocoon material anymore, namely, when 
numerical diffusion has lowered the entropy inside the co- 
coon to levels comparable with the ambient entropy. 

Numerical mixing is an irreversible process and hence 
creates entropy. On the other hand, 'mixing' across the con- 
tact discontinuity in real radio galaxies is not necessarily 



irreversible as can be seen from the following argument. 
Suppose that the ambient ISM/ICM material is not mixed 
on the microscopic level. Instead, assume that intact dense 
ISM/ICM filaments and clouds are simply swept up into 
the light cocoon and kept distinct due to being on different 
magnetic structures. The fact that this is a thermodynam- 
ically reversible process can be seen by noting that, given 
a sufficiently long time, buoyancy effects could separate the 
ISM/ICM gas from the cocoon plasma. Thus, the entropy 
structure of the simulated cocoon becomes increasingly less 
physical as time proceeds. Indeed, once numerical mixing 
has significantly affected a substantial part of the simulated 
cocoon (which occurs at f ~ 8), the entropy structure of the 
cocoon may be a poor reflection of reality. 



4-2.2 Definitions of energies 

To explore the energetics of our simulation in a quantitative 
manner, we compute the following quantities: 

(i) Total mass of cocoon. 
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rium with the ambient materiah After the Sedov phase, the 
cocoon source then enters the buoyant phase. 

At all times during the buoyant phase in our simu- 
lations, the rising buoyant plumes entrain and lift mate- 
rial (mainly in the trailing regions). This contributes to 
an increase in the gravitational potential energy of both 
the plume and ambient material. Concurrently, an approx- 
imately spherical pulse propagates radially outwards in the 
ISM/ICM atmosphere. This pulse is the remnant of the 
strong shock that bounded the supersonic cocoon at early 
times. As will be discussed below, this pulse mediates a gen- 
eral expansion (inflation) of the ISM/ICM atmosphere which 
results due to a heating of the core regions by the jet activity. 



Figure 5. Normalized volume distribution of material with given 
specific entropies (computed over the whole computational do- 
main). The solid line corresponds to very early times t = 0.02, 
the dotted line corresponds to t = 2 and the dash-dotted line 
corresponds to t = 4. The second broad peak in the distribution, 
which corresponds to the cocoon material, gradually moves to 
lower entropies as it is polluted by low entropy ISM/ICM gas. 



Mcocoon(t) = / pdV, (6) 

Jc 

(ii) Total internal energy of cocoon (ambient) material, 

-Bint,coc{amb)(i) = j- / P dV, (7) 

(iii) Total kinetic energy of cocoon (ambient) material 

-Bkin,coc{amb)(*) = ^ / pv^ dV, (8) 

^ JC{A) 

(iv) Total gravitational potential energy of cocoon (am- 
bient) material 

-Epot,coc{amb)(t) = - / dV, (9) 

where C is the region in which s > 10, and A is the region in 
which s < 10. We then reference all energies to their values 
at the initial time: 

AE ^ E{t)~ E{0). (10) 

4-2.3 The active and early-inactive phases 

The bottom panels of Fig. ^ shows how these energies vary 
with time (upto t — 10) for Run 4. During the active phase, 
the energy injected by the jets is transformed into internal 
and kinetic energy of both the cocoon and ambient mate- 
rial, as well as the potential energy of the ambient material. 
All of these energy forms increase linearly with time as the 
injected energy is shared roughly according to equipartition. 
After the source turns off, the kinetic energy of the cocoon 
rapidly decreases as the rapid motions associated with the 
jet cease. During the short Sedov phase that follows (be- 
tween t = 1 and t « 1.5), the internal energy of the cocoon 
also undergoes a rapid decrease (adding mainly to the in- 
ternal energy of the ambient medium) , as the overpressured 
cocoon expands adiabatically to achieve pressure equilib- 



4-2.4 The entropy structure at late times 

The radio galaxy influences the entropy of the ISM/ICM 
atmosphere by driving shocks into it. At very long times 
after the jet activity has ceased, the ICM/ISM will be in 
a relaxed hydrostatic equilibrium in which the entropy is 
increasing outwards (or else there will be convective turnover 
until such a state is achieved). A parcel of gas which passed 
through a very strong radio-galaxy induced shock (which 
may have initially been in the centralmost regions of the 
galaxy /cluster) will have a high final specific entropy and 
achieve its new equilibrium state some large distance from 
the cluster core. If the same parcel of gas suffers a weaker 
shock, the entropy increase will be less and it will settle at 
a lower point in the final ISM/ICM atmosphere. 

The density, temperature, and entropy plots in Fig. ^ 
show that, after the episode of radio galaxy activity, the 
ISM/ICM core is left on a higher adiabat than it was prior 
to the activity. This is easy to understand. Since the ini- 
tial radio galaxy expansion was supersonic, the material 
which initially had the lowest specific entropy (i.e. at the 
center of the initial ISM/ICM distribution) has all been 
shock heated. The final ISM/ICM core contains a mixture 
of (mildly) shocked gas, and unshocked ISM/ICM which has 
adiabatically flowed inwards from more distant regions (and 
thus possesses a higher specific entropy) to displace highly 
shocked material. The highly shocked material is expelled to 
the outer regions of the galaxy/cluster. For the specific case 
represented by our simulation, the excess entropy over the 
undisturbed cluster gas peaks at 20%. This translates into 
an excess entropy of about 14 keV cm'^ . 

Although this is a small change in the entropy of the 
cluster core, one might think that repeated radio galaxy ac- 
tivity would lead to a significant accumulation of excess en- 
tropy in the cluster core. This depends upon the assumed en- 
tropy profile for the ISM/ICM. Within the /3-model param- 
eterization used here, the entropy profile of the ISM/ICM is 
rather shallow. Thus, even if one could imagine completely 
shocking a large part of the central ISM/ICM, the entropy 
increase of the core seen in the final equilibrium would still 
be rather modest since material just outside of the strongly 
shocked region will simply flow inwards and take its place. 
If the initial entropy profile is much steeper (as it is in some 
strong cooling flow models) then the same radio galaxy could 
impact the entropy of the core rather more signiflcantly. 
Clearly, further work is needed on the contribution of ra- 
dio galaxies to the ISM/ICM entropy in which the entropy 
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Figure 6. A comparison of the initial and final states of the simulation. Top left panel: Total energy contained within radius r at two 
late times it = 19 and f = 20) minus the same quantity evaluated at the initial time (t = 0). The solid and dotted lines show this 
energy-difference evaluated for times t = 20 and t = 19, respectively, for regions along the jet axes id < 30° and > 150°). The dashed 
and dot-dashed lines show this energy-difference evaluated for times t = 20 and t = 19, respectively, for regions away from the jet axes 
(30° < < 150°). The fact that these last curves are higher at r = 30 shows that most of the energy injected by the radio galaxy 
resides in the expansion wave at large angles from the jet axis. Top right panel : Solid and dotted curves show the angle-averaged density 
profile at t = 20 and t = 19, respectively, ratioed against the initial density profile of the ambient material for regions away from the jet 
axis. Middle left panel : Solid and dotted curves show the angle-averaged temperature profile at t = 20 and t = 19, respectively, ratioed 
against the initial (isothermal) temperature profile of the ambient material for regions away from the jet axis. Middle right panel : Solid 
and dotted curves show the angle-averaged entropy profile at f = 20 and t = 19, resijectively, ratioed against the initial entropy profile 
of the ambient material. The entrojjy enhancement in the core of the ISM/ICM atmosphere is clear. The large fluctuations within the 
core radius seen in these last three panels are due to the high-entropy blobs of cocoon/plume material which are rising slowly under the 
action of buoyancy forces. The vertical dashed line in each plot denotes the location of the core radius. Bottom left panel : Total energy 
contained within radius r at late time (t = 20) minus the same quantity evaluated at the initial time (t = 0), with each of the three 
forms of energy — gravitational potential energy (solid line), kinetic energy (dashed line) and internal^ ergrgg jotted g^Q^j^Q^Q^Q^^^Q 
separately. Bottom right panel : Total energy contained within a given mass coordinate m at late t^e^^= §dj itimus^^^sanie? quanfi^y 
evaluated at the initial time {t = 0). 
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injection is considered in relation to particular detailed cool- 
ing flow models. 

Using 7? 05^ r data, Ponman, Cannon & Navarro (1999) 
and Lloyd-Davies, Ponman & Cannon (2000) have found ev- 
idence for an "entropy floor" in groups and clusters. This 
suggests from an observational point of view that these 
systems possess excess entropy (as compared to the en- 
tropy profile expected from simple scaling relations and 
detailed simulations; however, see alternative viewpoint of 
Bryan 2000). Lloyd-Davies, Ponman & Cannon (2000) claim 
an excess entropy of 70-140 keV cm'^ in a sample of 20 
clusters and groups. Interestingly, these authors also con- 
clude that this entropy injection occurred at low densities 
(1 — 3 X 10"* cm~^ ) and hence they suggest that it occurred 
prior to the collapse/formation of the cluster. In principle, 
microscopic mixing of the ICM and radio plasma inside of 
the cocoon/plumes might also lead to ICM entropy injec- 
tion at low density. However, our simulations suggest that 
buoyancy effects tend to transport any such mixed material 
completely out of the cluster center. Thus, if the determina- 
tion of the density at which the entropy injection occurred 
is robust, it does indeed seem likely that the bulk of the 
entropy injection occurred prior to cluster formation and is 
not due to more recent radio galaxy activity. 



4-2.5 Energetics at late times 

The energetics of the very late stage of Run 4 are reported 
in Fig. ^ The large scale pulse mentioned in §4.2.3 manifests 
itself in each of the quantities plotted in Fig. ^ at r ~ 20. The 
top-left panel shows the difference in enclosed energies as a 
function of radius between late times and the initial state 
t = 0. Comparing the values of this function on each side 
of the pulse may suggest to the casual reader that the wave 
contains most of the injected energy. This is a curious result 
since we have already noted appreciable entropy and internal 
energy increases associated with the irreversible heating of 
the ISM/ICM atmosphere. 

The bottom panels of Fig. ^ elucidate the true nature 
of this pulse. In the bottom-left panel, we note that the 
major contributor to the increase of enclosed energy comes 
from gravitational potential energy (rather than kinetic and 
internal energy as would be expected if this were a sound 
wave carrying away the injected energy). Secondly, display- 
ing the enclosed energy difference in mass coordinates (Fig. ^ 
bottom-right panel) shows that the mass contained within 
the inner edge of the wave-packet has acquired an apprecia- 
ble excess of energy during the evolution corresponding to 
at least half of the injected energy. 

Thus, it can be seen that the ISM/ICM atmosphere has 
been inflated (or puffed up) due to the heat deposited in the 
core of the atmosphere by the jet activity. The spherical 
pulse (travelling at the sound speed) drives mass to a higher 
location in the potential thereby mediating this inflation. 
That this inflation has occurred is obvious from the reduced 
core density in the final state (Fig. |^ top right) . Examining 
the enclosed energy in mass coordinates lets one see that 
this infiation actually accounts for a substantial fraction of 
the injected energy, implying that the injected energy has 
been heated the ISM/ICM with a high efficiency (r? ~ 0.5). 



5 DISCUSSION 

5.1 Cooling flows and radio galaxy heating 

As discussed above, we find that a large fraction (77 ~ 0.5) of 
the energy injected by the jet activity is thermalized in the 
core regions of the ISM/ICM. As mentioned in the Intro- 
duction, this raises the possibility that radio galaxy activity 
can have signiflcant effects on the energetics of the ISM/ICM 
and, therefore, act as negative feedback for cooling flow ac- 
tivity. 

In particular, let us address the case of cluster cool- 
ing flows. Peres et al. (1998) compared the radio and X-ray 
luminosity of 58 clusters in a flux limited ROSAT sample. 
They found strong radio sources in about 10% of the clus- 
ters within the sample. However, using the results of Peres 
et al. (1998) and summing over all members of the sample 
(not just those with radio sources), we find that the ra- 
dio luminosity of the sample (defined as the sum over all 
sources of uLv measured at 1.4 GHz) is about 1% of the 
sample bolometric X-ray luminosity from within the "cool- 
ing radius" . To discuss the relevance of radio galaxy heating 
on the cluster, we must relate the radio luminosity to the 
kinetic luminosity of the sources. 

Simple dynamical models of expanding radio lobes sug- 
gest that the maximum synchrotron luminosity of homoge- 
neous, space-filling, lobes is only a few percent of the kinetic 
luminosity (e.g., see Bicknell, Dopita & O'Dea 1997). De- 
partures from equipartition always act to decrease the syn- 
chrotron efficiency, and the efficiency further decreases as 
the source expands (Begelman 1996, 1999). The only way to 
increase this efficiency is to concentrate the radio emission 
into a small volume of the radio lobe (e.g., strong shocks 
within the lobe). In the absence of such inhomogenieties, 
it seems likely that the radio luminosity is at most a few 
percent of the kinetic luminosity. Therefore, using the Peres 
et al. (1998) results quoted above, the sample-averaged ki- 
netic luminosity is comparable to the X-ray luminosity from 
within the cooling radius. The results of our work suggest 
that a large fraction of this jet power can be thermalized in 
the ICM. 

Thus, from a purely energetic point of view, radio galax- 
ies can balance the radiative cooling associated with cooling 
flows. However, it is clear that the ICM of real galaxy clus- 
ters does not experience a steady-state situation in which 
jet-induced heating balances radiative cooling. As has al- 
ready been noted, most of the integrated radio luminosity 
in the Peres et al. (1998) sample is contributed by the power- 
ful radio galaxies that reside in only 10% of the clusters. The 
arguments of the previous paragraph would suggest that, in 
these clusters, the radio galaxy heating would vastly exceed 
the radiative cooling. Conversely, radiative cooling would 
certainly seem to dominate in the remaining 90% of the sys- 
tems (unless they possess powerful, but invisible, nuclear 
outflows; Binney 1999). Thus, it is clear that even if radio 
galaxy heating is relevant, the heating/cooling balance can 
be vastly different from cluster to cluster, or at different 
points of a given cluster's lifetime. 

It is beyond the scope of this paper to confront, in any 
direct manner, the plethora of new data that the Chandra 
X-ray Observatory and XMM-Newton are providing on the 
properties of the ICM and cooling flows. Here, we simply 
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note some of the challenges that future models of cooling- 
flows and embedded radio-galaxies must meet. 

Chandra finds spectroscopic evidence for cooling flows 
(albeit with a mass deposition rate smaller than previously 
estimated with ROSAT; Peres et al. 1998) even in clusters 
that contain energetic radio galaxies. For example, Fabian 
et al. (2000) find cool gas in the Perseus cluster, even though 
the central radio galaxy is clearly interacting with the ICM 
is a rather vigorous manner. Curiously, the cool gas seems 
to form a shell around the radio lobes and occupies precisely 
the location where we would expect shock/compressionally 
heated gas to reside. Furthermore, Chandra has observed 
sharp temperature discontinuities in clusters that do not ap- 
pear to host powerful radio galaxies. Firstly, a few clusters 
(e.g., Abell 2142, Markevich et al. 2000) have large scale 
temperature discontinuities which have become known as 
"cooling fronts" . These may bo discontinuities in the cooling 
flow associated with a group/cluster merger event. Secondly, 
Mazzotta et al. (2001) has found a bubble of hot gas in the 
cluster MKW3S. They associated this with some previous 
outburst of nuclear activity from the central galaxy. How- 
ever, if such a bubble had buoyantly risen from the centre 
of the cluster, we would expect KH instabilities to disrupt it 
to a larger extent than observed. One would need to invoke 
some mechanism, such as large scale magnetic fields, to sta- 
bilize the large scale KH instabilities and prevent disruption 
of the bubble. 

XMM-Newton has also produced a further mystery. By 
performing detailed emission line spectroscopy of several 
cooling fiow clusters (e.g., Abell 1835, Peterson et al. 2001; 
Abell 1795, Tamura et al. 2001) with the reflection grat- 
ing spectrometer (RGS) on XMM-Newton, one can probe 
the detailed temperature distribution of the ICM. Multi- 
ple temperatures were observed in concordance with cool- 
ing flow models, with components ranging from the virial 
temperature down to ~ 1 keV. However, no emission lines 
were observed from gas below ~ 1 keV whereas cooling mod- 
els would predict large amounts of emission from such gas. 
Fabian et al. (2001) discusses various solutions for this puz- 
zle and favours a model in which there are large abundance 
inhomogenieties affecting the predicting cooling rates. An- 
other solution is to postulate a heating mechanism which 
targets cool (~ 1 keV) gas. Examples of such a mechanism 
are the "percolation" model of Begelman (2001), or the "re- 
connection" model of Norman & Meiksin (1996). Heating 
models based upon simple hydrodynamic shocks (such as 
presented in this work or that of David et al. 2001) do 
not possess such temperature dependent targeting. In or- 
der to match the observations, the putative heating mech- 
anism must meet the further challenge of heating the cool 
gas rapidly up to the virial temperature. 

Finally, we note that neither Chandra nor XMM- 
Newton have yet found clear evidence for strong shocks asso- 
ciated with radio galaxies. Using the characterization of Pa- 
per I, all cluster radio sources observed so far appear to be in 
the weak-shock or the no-shock regime. Although more work 
is required before making a rigorous statement, there ap- 
pears to be a disagreement between the observed frequency 
of strong shocks and current hydrodynamic models of radio 
galaxy evolution. This might suggest that radio jets are less 
powerful than we believe (i.e. the synchrotron efficiency is 
higher than a few percent, implying that the radio emitting 



regions fill a small portion of the radio lobes). This problem 
was first pointed out by David et al. (2001). 

5.2 The limitations of this work 

The simulations that we present in this paper are just the 
first step in addressing the long-term evolution of radio 
galaxies. In order to make the problem tractable, we have 
made several simplifying assumptions. Here, we briefly dis- 
cuss how these assumptions might affect our results. 

Our biggest simpliflcation is that we have reduced the 
problem to two spatial dimensions by assuming strict ax- 
isymmetry about the jet axis. Given the inherent axisym- 
metry of the system, and the fact that observed radio co- 
coons are approximately axisymmetric, this is not an out- 
rageous assumption. However, there are several undesirable 
consequences of assuming axisynunctry. Firstly, we prevent 
the formation of non-axisymmetric KH and RT modes. All 
of the available free energy is then channeled into the ax- 
isymmetric and RT modes. The axisymmetric KH modes are 
most efficient at mixing (since they present the most effective 
barrier to the flow), and hence the axisymmetry assumption 
maximizes the amount of mixing between the ambient and 
cocoon material. Secondly, the assumed axisymmetry pre- 
vents us from including the effects of any realistic random 
ISM/ICM motions (i.e., "weather"). These ICM motions will 
buffet the buoyantly rising plumes and will become crucial 
in determining the morphology of very old sources. 

We have also chosen to ignore the presence of magnetic 
fields (i.e., we perform a hydrodynamic, rather than a mag- 
netohydrodynamic, calculation). Magnetic fields in the radio 
cocoon might well be in equipartition with the particle en- 
ergies and hence could be dynamically significant. Dynami- 
cally important fields that arc tangled in an isotropic man- 
ner on small scales behave like a relativistic gas and so can 
be handled within the context of a pure hydrodynamic sim- 
ulation. However, anisotropic tangled fields, or large scale 
magnetic fields, will introduce new features to the dynamics 
which require a full magnetohydrodynamic treatment. Sim- 
ilarly, the importance of neglecting the cocoon's relativistic 
equation of state can only be assessed once calculations are 
performed using relativistic hydrodynamic codes. 



6 CONCLUSIONS 

Radio galaxies are known to be dynamical systems that 
evolve on timescales of 10^ — 10* yr. Hence, we expect many 
galaxies and clusters of galaxies to be host to dead, or relic, 
radio sources. It is possible, and maybe even probable, that 
this past radio activity has influenced the energetics and 
thermodynamics of the hot, space-flUing medium in these 
hosts. Thus, a study of the environmental impact of dead 
radio galaxies is important and, given the high quality X- 
ray data now coming from the Chandra X-ray Observatory 
and XMM-Newton, very timely. 

We have used high resolution hydrodynamic simulations 
to investigate one particular scenario relevant to dying radio 
galaxies — i.e., a rather powerful radio- loud AGN situated 
at the center of a galaxy or cluster of galaxies whose activ- 
ity ceases abruptly. In particular, we have simulated back- 
to-back jets propagating in a /3-model galaxy/cluster atmo- 



© 0000 RAS, MNRAS 000, 000-000 



The hydrodynamics of dead radio galaxies 11 



sphere. We then shut down the jet activity and let the result- 
ing structure evolve. To make the problem more tractable, 
we assume axisymmetry about the jet axis (thereby reduc- 
ing the problem to two spatial dimensions) and neglect the 
action of magnetic fields. This is a significant extension of 
the work of Churazov et al. (2000) and Briiggen & Kaiser 
(2000) since we simulate the active phase of the source and 
then shut off the jets, rather than modelling the late stages 
of evolution by simply letting a static bubble evolve under 
the action of buoyancy. 

In the early lifetime of the simulated source (while the 
radio jets are still active) , the shocked jet material forms a 
cocoon which is bounded by a shell of shocked ICM. Even 
during the active phase, hydrodynamic instabilities start 
to shred this cocoon. However, only after the jet activity 
ceases do KH and RT instabilities destroy the integrity of 
the cocoon. Thereafter, the old cocoon material forms two 
plumes which rise in the cluster potential through the ac- 
tion of buoyancy forces. These plumes entrain a significant 
amount of cooler material from the ICM core and lift this 
material high up into the cluster atmosphere. At very late 
times, the galaxy /cluster atmosphere settles back into hy- 
drostatic equilibrium but with a core specific entropy that 
has been enhanced by ~ 20% over its initial value. This 
entropy enhancement is due to the shocking of the lowest 
entropy material by the strong shock which bounds the ra- 
dio cocoon during the early active phase. We find that a 
large fraction (~ 0.5) of the injected energy is thermalized 
in the ISM/ICM. Comparing late times to the initial state, 
we find that the ISM/ICM atmosphere has become inflated 
in order to maintain hydrostatic equilibrium after the ther- 
malization of the jets energy. Such a large thermalization 
efficiency raises the possibility that radio galaxies are im- 
portant in the overall energy budget of the ISM/ICM. 
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APPENDIX 

In this appendix, we present a brief derivation of the relativistic dispersion relation for the Kclvin-Hclmholtz instability in 
cylinderical geometry. We then use this dispersion relation to connect our simulations to the behavior of real sources. Similar 
relativistic KH analyses have been performed by many authors (Blandford & Pringle 1976; Ferrari, Trussoni & Zaninetti 1978; 
Ray 1981; Birkinshaw 1984; Zaninetti 1985). Wo repeat the derivation hero for completeness. 
The continuity and momentum equations of relativistic fluid dynamics for a 7-law gas are 

irMl = (11) 

ox" ax^ ox'' 

where p and p are the comoving density and pressure, respectively, w' are the components of the four velocity, and 7 is 
the polytropic index. Here, /i = p + 7p is the relativistic enthalpy. We wish to use these equations to study the relativistic 
Kelvin-Helmholtz instability. 

Consider first-order perturbations of these equations: 

p = po + Ap (14) 
p = po + Ap (15) 
u' = ul + Au\ (16) 

For simplicity, we shall drop the subscript '0' in the subsequent equations. Assume that the unperturbed quantities are 
time-independent and spatially uniform. The linearized fluid equations can then be cast in the form 

h A d^Ap _ d''Ap 

7P J dx''dx* dx^dxi 

For definiteness, we choose cylinderical polar coordinates and say that the unperturbed velocity fleld is along the 2:-direction 
with velocity V . This last equation can then be written as 

Now consider the following geometry which is appropriate for the analysis of the Kolvin-Hclrnholtz instability in radio 
galaxy cocoons. Consider a contact discontinuity situated at the surface r = R. For r < R, we have a fluid (denoted by 
subscript 1) which has an unperturbed velocity Vi in the ^-direction. For r > R, we have a fluid (denoted by subscript 2) 
which has an unperturbed velocity of V2. We search for perturbed modes of the form Ap = /(r) exp[i{kz + m(f> — ujt)], imposing 
a regularity condition at r = and only permitting outgoing waves in the region r > R. Matching the pressure perturbations 
at r = i? gives the following dispersion relation: 



V p,rifel(cj- Wl)-'= V"- V V.^h2(uj - kV2)\ (19) 

Kx3'm[KxR) K2H^r^^'{K2R) 

where 

Kl = vl(^- 1] (a; - kVif - k"" + a;^ (20) 
V71P / 

with a similar expression for K2. Without loss of generality, we can suppose that the fluid in the region r > i? is at rest 
(V2 = 0, Vi = V). It can be seen from this last expression with a little algebra that, in limit of small wavelengths and only 
mildly relativistic wave speeds, the general solution has the form 

a; = fc/(/ii//i2,X,7i,72) (21) 

where M. = V/cs,2, where Cs,i is the sound speed in the region r < R. 

Wc can now use this general solution to comment on the connection between our simulations and the hydrodynamics of 
real sources. Since most of the interesting hydrodynamics is associated with the development of the KH instability, we should 
endeavour to approximate the expected real-life KH growth rates within our simulations. 

The cocoons in real radio galaxies will, at least initially, be dominated by plasma with a relativistic equation of state 
(i.e. 71 = 4/3). On the other hand, our simulations are inherently non-relativistic (i.e., 71 = 5/3). While fully relativistic 
simulations are needed to address, in full, the effect of a relativistic equation of state, we can use the above dispersion relation 
to judge the effect on the KH growth rates of neglecting the relativistic nature of the cocoon plasma. Numerical solutions of 
eqn (19) show that the KH growth rates and wavelengths are little affected by changing from 71 = 4/3 to 71 = 5/3. Hence, 
neglecting the relativistic equation of state is not a serious flaw of our simulations. 

Noting that the Mach number of the cocoon backflow is determined by the hydrodynamics of the system (and is always 
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approximately unity), the remaining parameter relevant to the KH growth rate is the ratio /11//12. In real systems, the 
ambient thermal material will be non-relativistic (with sound speed 05,2) while the cocoon (especially near the jet head) will 
be dominated by relativistic plasma (hi « p). Since c ^ Cs,2, the ratio of relativistic enthalpies will be 

¥ = ^ (22) 

/l2 p2C^ 

(23) 



Plausible parameters give Cs_i ~ 3000 — 4000 kms ^ , giving hi/h2 ~ 1 — 2 x 10 Any pollution of the relativistic plasma in 
the cocoon (due to mixing) will reduce this ratio appreciably. For comparison, in our simulations we have 

= - (24) 

/l2 P2 



2 

Cb,1 



(25) 



and by measuring the sound speeds, we determine that we achieve h\/h2 ~ 5 — 10 x 10~ . We conclude that while it would 
be desirable to push this ratio lower by a factor of a few, our simulations are exploring the appropriate parts of parameter 
space. 
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